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The problem of studying the growth of the error is most important for the numerical solution of, 
differential equations. In this paper the Wilfs criterion is generalized to be applied for systems of 
differential equations. A general theorem is investigated and regions of stability have to be deter- 
mined. The use of an electronic computer is more essential for such regions to be characterized. These 
regions of stability have the property that, the error introduced at any stage tends to decay. The regions 
of stability for particular numerical methods are explicitly determined. 
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1. Introduction 

The system of m ordinary differential equations can be written in the vector form as follows: 



dx 



with the initial value Y(xo) — Fo, 
where 




F(x;Y)A 



/fi(x; y u y 2 , • . ., ym) < 
f 2 (x; y u y 2 , . . ., ym) 

K fm(x; y u y 2 , . . ., y m ) * 



Let the vectors: 

Y(x) be the exact solution of the system (1), 

Y n be the theoretical (approximate) solution at x = %m 
Y n be the practical (computed) solution at x = x n . 

We define the error e M of the solution at the point x n by: 



€n = Y n -Y(x n ). (2) 



*An invited paper. 
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Hence we get, 



* 



e' n = Y' n -Y'{x n ) 

=F(x n ; Y n )-F(x n ; Y n — e n ). 

Assume that: 

(a) F{x n \ Y) is a continuous function of Y for Y contained in the closed intervals whose end- 
points are Y(x n ) and Y n , 

(b) the matrix ( — -r-7 J; (£,7=1, 2, . . ., m) exists for Fin the open interval whose end- 
points are Y(x n ) and Y n . 

Therefore, by the mean value theorem we get: 

€n= \d) en = Jn€n > ( l "^~ 1 ' 2 > • • •' OT )' ( 3 ) 

where the elements of the Jacobian matrix J n = (t~M are taken at suitable places 

Wj/n 

(x ny Fn-diag (6i) -€ n ),0<ft<l, £=1, 2, . . ., w. 

THEOREM: Stability criterion for numerical integration methods for the solution of systems 
of differential equations. 

The integration formula 



Y n+1 = J a„Y n+ „+h j b,Y; + , (4) 

/or £/ie solution of m ordinary differential equations Y' = F(x; Y) is stai/e i/"a/id o/i/y £/, £/ie Hermi- 
tian forms with the coefficients 






, minors) - _ 

A(rs)= 2^ l^q+1+l-r, i^q+1+l-s, i — C r _i 5 jC s _ l3 ij, 



(5) 



/or all i = 1, 2, . . ., m are positive definite 

where Cg, 4 = K-jb^-q — ag_ q , 8 = 0, 1, . . ,,q+l (5) 

with K i = — hp i , h the step of integration, ai=— 1, 

pi = the eigenvalues of the Jacobian matrix ( — L 

1 = an extrapolation formula 

if (4) 
7^ an interpolation formula. 

Proof: The formula (4) gives the theoretical approximating vector y n +i for the exact value 
F(x w + A). The computed (practical) solution at x = x n + i is determined as follows: 
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y„ +l = J a„Y l ,,„ + h j 6,.F;, +1 . + «„ +i , (7) 

where /?»+i the round off error. 

The corresponding exact solution is 

Y(x* + h)= X ajr(xn + vh)+h £ W(*ii + i*) + 7V+i, (8) 

where the truncation error 7*+i = lto+ih k+1 Y[*+p + 0(h k+2 ) 9 k ^ 1, yLC A + i # 0. Subtracting (8) from (7) 
we have 

F„ +l -F(*„ + /i) = 2 o„(F„ + „-Ffe + ^)) + /i 2 6»(I^i + ,-I"(*i + »A))+il, (9) 

where /?„ + i — 7a-+i = /? is generally a nonzero vector. 

Substituting (2) and (3) in (9) we get the difference equation for the error 

1 

€«+!= 2 (lv€n+v + h 2 bJ n+ r€n+i> + R, (10) 

where Jn+y = 1 T~~ J n + f . 

In studying the local growth of the error e n it is reasonable to assume that R and the matrices 
dfi\ in a small interval (x- Q , %\) are both constants and they may be written in the form R and 

9yJ 

J respectively. In practice they vary slowly from step to step, otherwise the step size of inte- 
gration is too large, Hamming [2J. 1 Under these assumptions we get from (10) the nonhomogenous 
system of linear difference equations with constant real coefficients in the form: 

1 

€*+l = 2 Clven+p + hJ ^ bvCn+p + R. (11) 

v=-q v=-q 

To solve this system, we introduce a new vector variable e in place of the variable e by means of a 
nonsingular linear transformation: 

e=T l e, where det T ^ 0. (12) 

Therefore the system (11) is transformed into: 

i a 

€n+i= 2 a v € n +p+hK ^ b v € n +v+R, (13) 

v=-q v=-q 

where 

T i R = R, K=T l JT. 

In particular the transformation T may be chosen so as to put the matrix T^^T into the classical 
canonical form, Zurmiihl [3J. 



1 Figures in brackets indicate the literature references at the end of this paper. 
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The corresponding homogeneous system of (13) has the form 



€n+i= 2 dvtn+v + hK 2 ^*n+ 



(14) 



Case I: If the eigenvalues pi of the matrix J are all distinct then the matrix /£ has these eigen- 
values in the principal diagonal and zeros elsewhere. In this case the system (14) has the form: 



(Kj5l-f l)€i,n+l= ^ ( a v~ K M*i,n+v, 
1=1, 2, . . ., 771, 

where 

Ki = — hpi. 

Put €i, y = X£, Norlund [4], we get from (15) the characteristic equations: 



1 = 1, 2, . . ., m. 



(15) 



(16) 



Case II: If the eigenvalues p, of the matrix J are not all distinct, then the matrix J may be 
transformed to: 



K = 



& 



K t 



.0 



Ki ei 



K,, 



K u 



where each submatrix Ku is of the form: 



Kn = 



Pi 1 

Pi 1 

1 







P> J 



with pi in each position in the leading diagonal, unity in each position in the diagonal immediately 
above it, and with zeros elsewhere. In this case of nondistinct eigenvalues, we have for every 
submatrix of order y a system of equations in the form: 
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6i, n+i 



€2, n-f 



€>, n+\ 



2 a„€i, //tr 

() 



V a^6y. n + i> 



+ /l 



Pa 



Pa 



1 
pa 



2 ^1, /H 



2 6A.« 



2 6^7, n + 



That means 



where 



v=-q v=-q 

* fl i<7 



£ = 7. 



£=1,2, . , 
The characteristic equations system can be got as follows: 

t'=() 

i=l,2 f . . ., 7. 

According to (16) and (18) the characteristic equations system has the form: 



Q 

1 



(K l -6, + l)Af+ 1 + 2] (Kib-u-a- v )kf- v = 



£=1,2, . . ., m, 



(17) 



(18) 



(19) 



independent of whether the eigenvalues of the matrix J are distinct or not. According to (6), the 
equations* system (19) has the form: 



2+i 



^here 



2a,iXr=o, 1=1,2, . . ., m, 



C„. i is generally a complex value. 



(20) 



Due to Wilfs definition [1]: The numerical integration method (4) is said to be stable if the 
roots of (20) for all i= 1, 2, . . ., m are inside the unit circle in the complex plane, so that the 
error introduced at any stage tends to decay than build up. 

In this connection we state a theorem of Schur [5]: 

The zero points of the polynomial: 

f(X)-=C + CiX+ . . . +c q+1 x q+1 
95 



are all inside the unit circle, if and only if the Hermitian form: 

H = 2 |C q +iXi + C q \ 1+ i + . . . +C 1+ i\ q | 2 |C Ai-f-CiA 1+ i+ . . . -hCq-^J 2 } 

1=0 

is positive definite.' 2 

The quadratic form H can be formed as follows: 

ff=iff; (ic<*.-ri»-ici»)iwi» 

1=0 lr=() 

+ 2R e 2j 2j (C Q+ i-jC q+ \- k — CjCk) X/+jX/ +A - [ 

j=0 k=j+l J 

=i{i(iQ + < + .-. s i 2 -ic s -,i 2 )iA S |* 

/=o U=z 

Q Q — _ _ ] 

+ 2/i e ^ 2^ (C^ + / + i_ r C 9+ /+i_ 6 . — C r -iC s -i)k r k s > 

r=l s=r+\ J 

=i w 2 i (ic, + , +1 - s i 2 -ic s ^i^ 

+ 2Rei J fAr\.2 (C, + , +1 _ r C, +/+1 _,-C r _,C,-,))- 

'•=0 s=r+l <■ ] =0 J 



With 



min (r, s) 

d rs = ^ {C7+/ + i- r C<7 + /+i-.«> — Cr-lCs-l}, r, S = 0, 1, . . . , <? 

f=0 



we have directly the form: 

Q - 

H= 2^ A rs \ r \s. 

r, s=o 

For the system of differential equations (1) we have the Hermitian forms 
#<*)= £ 4^ r X s , £=1, 2, . . ., m. 

r, s=0 

Therefore ^^ has the form (5) and accordingly the integration formula (4) is stable, if the Hermi- 
tian forms H^ are positive definite for all i, so that theorem (1) is proved. 

SPECIAL CASE: For m=l, we have one single differential equation in the form y' =f(x; y). 
The Jacobian matrix reduces to one single element df/dy and k = — hdf/dy is a real value. The 
stability criterion reduces to Wilf Stability Criterion [1]. 

2. Simple Examples 

Example 1: 

Y n ^ = Y n J rhY , tn T2=%h 2 Y'(x n ) +0(/i 3 ), (Euler-Cauchy method). (21) 



2 C means the conjugate complex value of 
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With g = 0, ao= 1, ai=— 1 and b {) = 1 we have from (5) and (6): 

^=|C,, < |*-|Co, i | 2 =l-|Ki-ll 2 . 

Thus A$ > 0, 1 = 1, 2, . . ., m, if and only if hp { is inside the circle \hpi + 1| = 1 (fig. 1) and so 
the integration method (21) is stable. 



f him (p.) 




Example 2: 



Figure 1 



1 4 2 



(22) 



r 2 =|& 2 r(*n)+0(A 3 ). 



14 2 

With g=l, c-i= — r, ao = T, ai=— 1 and &o = «, we get from (5) and (6): 



The determinants: 



4B=|G.,|Mc .,|»«4? f 



^-CuCrCuC^A 



1 2 

Co, i=T, Ci,i=— (k/ — 2) , C2, /= 1. 



n(0= j(/) = - 

^1 ^00 o' 



M°=(^!) 2 -W;, ) I 2 = ^{16-|3k ? --k / -4| 2 } 



With **=- (cu+ V-l /3/) it follows: 



M = -^{(a/ + 2) 2 + 4/3?-4}>0 

(a -h2) 2 
if and only if hpi is inside the ellipse -^ J — — - — h/3f = 1 (fig. 2) and so the integration method (22) 

4 

is stable. 
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i hIm(/>.) 

h/>.- PLANE 

-hRe(/o.) 



Figure 2 



Example 3: 



Y n +i = Yn + \ (Y' n + Y' n+l ) 9 Tz=-iY"(x n ) + 0(h*) 



(23) 



(Trapezoidal formula). 
With q = 0, Oo=l, fli=— 1 and 6o = 6i = ^ we get from (5) and (6): 

^ ( 8o=iCi,f| 2 — |Co,i| 2 =|§K t — 1| 2 — |i#Ci-hl| 2 = 2Re (k«). 

We have A$ > if and only if Re (p>) < 0, for h > 0, and so the integration method (23) is stable. 



Example 4: 



3. More Difficult Examples 



Y n+l =^Y n - i ^Y n ^h(2Y^Y; i+l ),T 4 = -^h^(x n )+0(h^ (24) 

14 4 2 

With 9=1, a_i=-, ao == F> ai= — 1> &o = 7 and b\=- we get from (5) and (6): 



45Hc,.i|»-|Co,i|»=><fl, 

C ,, = -i C 1(i =| (*-l), C,.,=| (2k, + 5). 



We have 



25 



^°=^I8 = ^{|2ici + 5|"-1}>0, 



if and only if ftp, is outside the circle |Ap/-2, 5| = 0, 5 (fig. 3), this means, the circle is the region 

of instability. 

Further: 

W={MY-\AW 

= ^{|2k,. + 5| 2 -1} 2 -^|(2k,- + 5)(k,--1)+k,--1| 2 



= z4{{|«'l 2 + 5Re (iti)+6} , -|2|j«i|* + 5it|-K|-6| t }. 

ozo 
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h hlmC^o.) 



hp.-PLANE 




|h,3.-2.5|=0.5 



Figure 3 

With Kt=- («/ + V^T Pi) we get &f >0 if 

{(a / -2)(a / -3)+/3?} 2 >4{{(a,-3)(a / +l)+/32}2 + 9/32}. 

The inequalities D$ > 0, i>=1, 2, are satisfied simultaneously, i.e., the integration formula (24) 
is stable if hp\ is inside the region R (fig. 4). 3 




Figure 4 



Example 5: 



Y M = Yn+±(-Y' m -t + BYl + Sr M ) 9 T 4 = -^Y"(x n )+0(h*).< 
With g=l, Oo=l, ai=— 1, 6_i=~T2, 6o=o and ^i=y^ we get from (5) and (6): 



(25) 



4flHc* f fMc ,il f 



= ^{|K/| 2 + 5Re(K»)+6}=^l, 

-4^1 = ^2, |Ci, i — Co, |Ci, i 



^(4|k/| 2 + 7*,-5ic,-12)=y«8, 



C*o,i i o K/, 



3 The regions of stability (figs. 4, 5, and 6) are determined by means of the electronic digital computer Z 22. 

4 This formula is given by Southard and Yowell [6]. 
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C u =j(2Ki-3). 



We have 



G,,= ^(5k, + 12). 



1 



W=6-Mkl 2 + 5ReM+6} 
=4{|k, + 2,5| 2 -0,25}>0 



if hp, is outside the circle \hpi~ 2,5|= 0,5 (fig. 3). 
Further: 

m={AW-\MY 

= ^{4{|K ( | 2 + 5Re (kY+W-WkV + Iki-SIc-UW. 
With Ki =- (aj + V^T /ft) we get D ( 2 » > 0, if 

{(ai-2) (a i -3)+ ) 8f} 2 >{(a i -2) (2a; + 3) +2/3?} 2 + 36/3f. 

The inequalities D^ > 0, i>=1,2 are satisfied simultaneously, i.e., the integration formula (25) is 
stable, if hp t is inside the region R (fig. 5) 



hlm(/»j) 




Figure 5 



Example 6: 

F„ +1 =-i Yn-2+l Y n - 1 +Y n +-^ (-4 r;_,+ii y;+5 k,; +1 ), ^=^ ^ n(x„) +o(*«) (26) 



(Hamming [2]). 
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With q = 2, a- 2 = — - , a-,=-, a = 1, «i=— 1, 

8 . 22 . 10 , 

b-, = — —, <>« = -^=, o,= — we get from (5): 

Co,i = g , Ci,j=— n^ (8k, + 3) , 

C 2j! = ^ (22k, -27), Cs..=^ (10k, + 27). 



The elements of the Hermitian Matrix follow 

A®=C 3 ,i C 3 ,i-Co,i Co,i=A% 

A < il=c z , i c u -Co, i c u =M, 

^11 ^3,1 ^3,1 + ^2,1 C2,i — Ci,i Ci,i —Co,} Coj. 

We have 

W=^=|C 3> i| 2 -|Co,i| 2 

1 



(27) J 



{|10/<i + 27| 2 -9} >0, if Zip, is outside the circle |Ap,--2,7| =0,3, (fig. 6). 



i\ hIm(/>.) 



h /> j - PLANE 



0.6 1.2 1.8 2.4V2.7//3 




hRe(^) 



|h/>j -2.7| = 0.3 



Figure 6 
Further 

^) = ^^-Mi/)| 2 = {|C3, 1 | 2 -|Co, i | 2 }{|C 3 ,v| 2 +|C 2 , j | 2 

- ICup-ICwl*} - \C 3 ,t C 2>i -Co,i C ljl | 2 = {|C 3> i| 2 - |C 0jj | 2 } 2 - \C U Cu-Cot c u \ 2 

= 7^fyT [25{5(af + #)-27 <* 1 + 36} 2 -4{{10(a? + /3?) -39 a,} 2 + 225/3?}]. 
From this it can be proved, that D.9' > for a,- < 0. Finally we get: 



A° = 



A® A<& M 

Aii m a\}} 

A\ii A\j} A® 



=4'K(^)*-Mlgl 2 }-2^l4?l 2 +2Re (4g C4</>) 2 ). 
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263-048 O - 67 - 



Put 



AfMAW-^m-iUrJ^,^), 



32 
'(27)« 



32 
(27)' 



-^^•-■^/.(oj.ft), 



— 32 

2Re(^(^)=^/,( aj ,ft), 



yhere 



/i(ai,/8i) = 4[65(a| + «)+87(yi + 180][25{5(af+/9f) 

-27a, + 36} 2 -4{(10(af + j 8?)-39a0 2 + 225/3?}], 
/«(ai,/ai)=-20[5(af+^)-27ai + 36][{55(^+j9f) 

- 87a* -180} 2 + 44100/3*], 
/ 3 (a i9 i3 i )=-8[(10(af + ^)-39a i ){(55(^ + )8f)-87ai 

-180) 2 -44100 ) 8f} + 6300 i 8f(55(a 2 + ^)-87a— 180)]. 
From these we get D\p > 0, if 

£(<*, fr) +/,(«■, fit) +/ 8 (ai f fid > 0. 

The inequalities D$ > 0, j/=l, 2, 3 are satisfied simultaneously, i.e., the integration formula is 
stable if hp x is inside the region R (fig. 7). 




Figure 7 
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